A ring of BEC pools as a trap for persistent flow 
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Mott insulator - superfluid transition in a periodic lattice of Josephson junctions can be driven by 
tunneling rate increase. Resulting winding numbers W of the condensate wavefunction decrease with 
increasing quench time in accord with the Kibble-Zurek mechanism (KZM). However, in very slow 
quenches Bose-Hubbard dynamics rearranges wavefunction phase so that its random walk cools, 
decreases and eventually the wavefunction becomes too cold to overcome potential barriers 
separating different W . Thus, in contrast with KZM, in very slow quenches is set by random 
walk with "critical" step size, independently of tq. As our study requires use of the truncated 
Wigner approximation (TWA) over relatively long time intervals, we investigate validity of TWA 
by comparing its predictions with exact calculations for suitably small quantum systems. 

PACS numbers: 03.75.Kk. 03.75.Lm 
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I. INTRODUCTION 

Gapless quantum critical points are a serious obsta- 
cle for quantum simulation with ultracold atomic gases 
or ion traps, where one would like to prepare a simple 
ground state of a simple initial Hamiltonian and then 
drive the Hamiltonian adiabatically to an interesting final 
ground state. This general observation has been recently 
substantiated by a more quantitative theory [l| (see [H 
for reviews) which is a quantum generalization of the 
classical Kibble-Zurek mechanism (KZM) Q. The the- 
ory predicts that density of excitations (or excitation en- 
ergy) decays with (usually a fractional) power of quench 
rate 1/tq. Experiments related to the quantum theory, 
although in the sudden quench limit tq — >■ 0, were made 
in Rcfs. 4, 5]. In one of them Q a ring of iV = 3 iso- 
lated Bose-Einstein condensates (EEC's) was connected, 
see Fig. [TJ Random initial phases of EEC's often result 
in a non-zero vorticity trapped in the final state. A sim- 
ilar ring with a bigger N can be prepared by painting_^a 
time-dependent potential, as is demonstrated in Ref. [6| . 

We use the truncated Wigner method in the Eose- 
Hubbard model to simulate gradually connecting a ring of 
condensates. For sufficiently slow quenches we find that 
typical winding numbers trapped in the connected ring 
do not depend on how slowly EECs are connected. The 
process responsible for the final trapped winding numbers 
can be thought of as frustration - trapping of the phase 
of the cooling, equilibrated condensate wavefunction in 
one of the minima corresponding to integer W. 

In contrast to KZM (where winding number W de- 
pends on the time t when the critical slowing down ceased 
to suppress phase-ordering dynamics [13), for sufficiently 
slow quenches the condensate phase evolves ergodically, 
exploring different potential minima. Therefore, instan- 
taneous half-width of a Gaussian distribution of winding 
numbers is determined by the step size of the random 
walk and proportional to the square root of the number 
of sites - its length. Trapping of this equilibrium winding 
occurs as the size of random walk steps - its "tempera- 



ture" - becomes too small to get over potential barriers. 
Final W are set by the critical step size CTc - the least 
step size that allows for jumps - and yield k, a^N . 

Thus, sufficiently slow quenches lose memory of phase 
dispersion and of W at t, which in faster quenches was 
shown to lead to consistent with KZM 12]. In the 
units (set by the dimensionless Hamiltonian immediately 
below) this (wide!) border turns out to be near tq ^ 10^. 

This paper starts with a brief dicussion of the Eose- 
Hubbard model (the subject of our study) in Section [H] 
and of the truncated Wigner approximation (our prin- 
cipal tool) in Section Hill Eefore presenting our results, 
we digress to explore the validity of TWA by comparing 
its predictions to the exact calculations of small quan- 
tum systems. This is done in the extensive Appendix. 
After we have introduced the model and validated the 
method we discuss results of our simulations and pro- 
vide an analytic understanding and a heuristic picture of 
the winding number trapping in very slow quenches in 
Sections iHlVini We conclude in Section HXl 



II. 



BOSE-HUBBARD MODEL 



The model describes spinless cold bosonic atoms in an 
optical lattice [Hi- In dimensionless variables, its 
Hamiltonian reads 



H = -J 
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with periodic boundary conditions. Here is a number 
of lattice sites and n is average number of atoms per site. 
For an integer n, there is a quantum phase transition 
from the Mott insulator to superfiuid at Jcr — 0. 
We drive the system by a linear quench 



J{t) = t/TQ , 

starting in the Mott ground state at J = 0, 
\n,n,n,...,n) , 



(2) 



(3) 



2 



J 



> J 



J 



BEC J tunneling 



where 9s = arg(0s) and each phase step 9s+i — 0s between 
nearest neighbor sites is brought to the interval (— tt, tt) 
modulo 27r. 

Figure [2] shows average square of the winding number 
W^, equal to its variance, at a final J as a function of 
the quench time tq. One can distinguish three different 
regimes of tq: 

• For small tq there is not enough time for phases 
at different sites to become correlated, the final (ps 
remain close to the initial fields and the random 
initial phases result in large winding numbers of 
variance = N/12; 



FIG. 1. An initial ring of isolated Bose-Einstein conden- 
sates, each with a definite large number of particles n and 
indefinite phase, becomes phase-correlated by slowly turning 
on the tunneling between the condensates, see Refs. ifii]. 



with the same large number of particles n at each site, 
and ending in the Josephson regime 

J < 1 (4) 

where the interactions still dominate over the hopping. 



III. TRUNCATED WIGNER 
APPROXIMATION (TWA) 

For large density n we replace annihilation operators 
Og by a complex field (f)s, Os ~ ^/rift's, normalized as 
^s=i I'^sP = N and evolving with the Gross-Pitaevskii 
equation (GPE) 



. d(j)s 



20, +(/.,_!) 1)0, . (5) 



Quantum expectation values are estimated by averages 
over stochastic realizations of 4>s. Each realization has 
different random initial conditions coming from a Wigner 
distribution of the initial state ([3|) : 



,(0) 



(6) 



with independent random phases 6',(0). 

TWA 0, II] is a semi-classical approximation accurate 
for sufficiently short quench times tq. In the Appendix 
we extract the largest tq where it is still applicable from 
simulations in small systems of a few lattice sites, and in 
Eq. (|23)) below we give an instanton-based estimate for 
large system sizes. Both estimates predict the largest tq 
to grow with the density n. 

Here we focus on the integer winding number 



W 
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• For longer tq there is more time to correlate nearest 
neighbor phases and, in accordance with the KZM 



[l2| - 14 j , W"^ decays with increasing tq . One might 
expect this KZM decay to extrapolate to tq oo. 

Contrary to this natural expectation, we find that 
for even longer tq there is a crossover quench time 
Tq where the variance saturates at a finite value 
which is the same for all J. Quite surprisingly, 
above Tq the winding number does not depend on 
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FIG. 2. Results from simulations using truncated Wigner 
method (|5I6[) averaged over 2 10* realizations. The plot shows 
variance of the winding number on N = 512 sites as a 
function of tq. For each J there is a threshold quench time 
tq{J) above which the variance saturates. The saturated 
variance does not depend on J. 



IV. JOSEPHSON EQUATIONS 

This unexpected result can be more readily explained 
in terms of the Josephson equations. We can always write 
l + /,)e*^= in Eq. © with real /s,6',. In the 
Josephson regime we have \fs\ ^ 1, see Fig. [3J3. After 
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FIG. 3. Results from simulations using truncated Wigner 
method (|5I6|) averaged over 2 10* realizations: In A, the satu- 
rated variance at J = 0.01 and tq = 52428.8 as a function 
of lattice size N. The solid line is a linear fit = 0.0079A'' 
for > 16. At A'^ = 8 the variance is below the linear fit 
and at = 4 (not shown) it is zero. In B, histograms of the 
modulus \4>s\ at J = 0.01 for tq = 52428.8 and = 512. 
In the Josephson regime J <S 1 fluctuations of the modulus 
around \(j}s\ — 1 are small. In C, a generic histogram of phase 
steps 61^-61^+1, here at J = 0.01, tq = 52428.8 and A^ = 512. 
The solid line is a Gaussian fit. 



elimination of fs we obtain Josephson equations 



2J(t) [sin {0, 



+1 



-i)] • (8) 



In case of constant J and more than 2 sites these equa- 
tions are chaotic. 

It is convenient to eliminate tq by introducing a 
rescaled time variable, 

U = t Tq'/' = Jit) t'J' , (9) 

when we obtain dimensionless chaotic equations 

d'^e.. 



du^ 



2u [sin {0,+i - 6s) - sin (0, - (10) 



with random initial phases 0^(0) and, for large n, vanish- 
ing initial velocities ^(0) = 0. 



V. THERMALIZATION 

The (dimensionless) nonlinear system (jlOp approxi- 
mately thermalizes after rescaled time u ~ 1 in the 
sense that averages of local observables can be obtained 
from a Boltzmann distribution. Since the Hamiltonian is 
time-dependent (jll2p . the temperature depends on time 
and the thermal distribution is only approximate. Nev- 
ertheless, in first approximation the evolution after u 
can be considered an adiabatic process with the state 



of the system following closely to the instantaneous state 
of thermal equilibrium. 
After u the variance 



(11) 



of phase steps Og+i — 0s is shrinking, see Fig. |4]4,, and en- 
ergy of the Josephson system becomes approximately 
quadratic 



s—l \_ ^ 




+ 2J{t) [1-cos {Os+i-es)] 



+ Jit) i9s+i - OsY 



(12) 



Consequently, the Boltzmann distribution of phase steps 
Os+i — Os becomes a Gaussian, see Fig. [3p, of zero mean 
and a variance related to the temperature 



T 



2JCT^ 



(13) 



of the quadratic system ([T^ . Here the Boltzmann con- 
stant fcs = 1. 

Due to cquipartition, thermal average of the quadratic 
energy (fT2l) is 



iE) 



T N 



(14) 



On the other hand, the time-dependent J in p2|) makes 
the thermal energy time-dependent 



dt^ ' ^ dt J 



(15) 



where (-Ekin) = ^TN is kinetic (hopping) energy, i.e., the 
thermal average of the last term in Eq. ((T2|) . A combi- 
nation of Eqs. (|14|) and ([TSl) gives a simple differential 
equation ^ ^-j or, equivalently. 



(16) 



This is the time-dependent temperature in the adiabatic 
process after u. 

A missing multiplicative integration constant in Eq. 
()16|) can be determined by an approximate initial condi- 
tion that (£') ~ JN at u. Here J ~ Tq corresponds 
to M ~ 1 and the initial condition assumes that at u 
the phases are (almost) as random as in the initial Mott 
state. With this initial condition at ii the temperature 
in the following adiabatic process becomes 



T 



J1/2J1/2 



^1/3 



(17) 



With Eq. (fT3)) this equation translates to 



-1/3 



-1/2 



(18) 



after u ~ 1. This scaling, in an equivalent form a 
is confirmed by numerical results in Fig. 2] A. 



u 

FIG. 4. In panel A, estimated dispersion u of a phase step 
between nearest neighbor lattice sites as a function of u on a 
lattice of A'' = 2048 sites for three different tq . The three plots 
collapse confirming that the rescaled time u is the relevant 
time variable. In the range 2 < u < 120 the collapsed log-log 
plots are linear with a slope —0.28 implying a ~ In 
panel B, corresponding evolutions of the winding number. W 
is trapped after a shrinks below CTc at Uc — 20. 



VI. ERGODICITY BREAKING 
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FIG. 5. Integrated phase 9s = X)j=i ((^i+i - ^j)\e(-^,n) 
as a function of the site number s and the rescaled time u. 
Notice that W = On/^tt. At u ~ 33.07 the winding number 
jumps by —1 at the link between sites 355 and 356, i.e., the 
phase step (0356 — 0355)\^(^_^ jumps from +7r to — tt across 
the saddle separating different W. 



The accuracy of the quadratic approximation in (jl2p 
suggests that with shrinking cr^ the system crosses over 
from ergodic to increasingly regular behaviour. The most 
striking manifestation of the increasing regularity is er- 
godicity breaking between different integer values of the 
winding number W, see Fig. |3|3. In energy landscape 
picture, valleys with different integer W are separated 
by an unstable saddle-point "phase slip" solution. In the 
Josephson regime J ^ 1 the phase slip is localized on 
a single link with a phase step 6s+i ~ Og = ±7r. The 
frequent initial jumps of W shown in Fig. |3}3 all pass 
through the localized phase slip, see an example in Fig. 

El 

According to the LAMH theory [H HI], when 
|VF/A^| ^ 1 the frequency of the winding number jumps 
is cx e~^'^"' , where the 4 J is energy (fT2)) of the localized 
phase slip. Since the temperature is /3^^ ~ 2Ja^ in Eq. 
(fT3)) . this activation coefficient is oc e'"^!" and the inte- 
ger winding number freezes out when a falls below 



1 



(19) 



From the data in Fig. |3}3 this happens at time Uc — 
20 — 50. This is when ergodicity between different Ws 
breaks down and the winding number (distribution) gets 
stuck. 



VII. 



TRAPPED WINDING NUMBER 



Below (Tc the phase step dispersion a keeps shrinking, 
so the random walk Qsiu) keeps smoothing, but its tem- 



perature is not enough to induce jumps unwinding its net 
winding number. Thus the frozen winding number is a 
remnant of a random walk with a = Or and its variance 



N at 



(20) 



Here Nd^^ is average distance-squared "random walked" 
by the phase around the ring. The estimate (pn)) agrees 
with the linear mW^ ^ 0.00797V in Fig. when 



0.56 



(21) 



For large N the linear scaling (I^Hl) gives stronger winding 
than typical winding originating from quantum fluctua- 
tions in the ground state that is only logarithmic in N. 

The winding number assumes a fixed value at a 
rescaled time Uc when a falls below Cc- For a given final 
J, as in Fig. [2l Uc translates to a quench time 



(22) 



which is independent of the lattice size N and whose 
dependence on J is consistent with Fig. [51 When tq ^ 
Tq the variance saturates at the finite value in Eq. 
(f20| . W settles down at Jc ~ UcT. 
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Q 



which is in the 



Josephson regime, Jc ^ 1, for slow enough quenches with 

^-^ 3/2 
TQ > Uc' . 



On the other hand, for tq <C u, 



3/2 



10'^, the familiar 



KZM scaling 



-1/3 



was observed [12| in quenches 



5 



that take the system beyond the Josephson regime into 
J 3> 1 territory. 

After the thermal freeze-out at Jc the winding number 
could still change by quantum tunneling. Its rate V cx 
with a numerical constant a ~ 1 can be obtained 



from instanton calculations 
TiyfTc ^ 1 or, equivalently, 



jig . It is negligible when 



TQ <^ V? . (23) 

With, say, n = 100 particles per site this is a generous 
upper estimate for the range of tq when the semiclassi- 
cal TWA is applicable. The same estimate p3| is also 
obtained from the condition that the winding number 
freezes out at a value of tunneling rate far above the Mott 
transition: Jc 3> n~^. Indeed, near the Mott transition, 
where the discrete nature of site occupation numbers is 
essential, the semiclassical TWA is not applicable. 



VIII. PHONONS 

After the winding number freezes out, it is convenient 
to think about a smooth persistent flow with phononic 
fluctuations on top of it. Pushing tq beyond Tq (or the 
rescaled time u beyond the freezing time Uc) makes no 
difference for the frozen winding number W , but it does 
make a difference for the phonons. Soon after the freeze- 
out most energy is deposited in the phononic fluctuations, 
but as a keeps shrinking below CTc, the field tends to a 
smooth persistent flow 



(24) 



where W is the winding number frozen near Uc and the 
circulation has constant phase steps Qs+\ — Os = 2t:W/N. 

The regular behaviour in this regime can be described 
by small phase fluctuations Vs (phonons) around this 
smooth background Og = 2ttWs/N satisfying a linearized 
version of the Josephson equations ([TU|. Indeed, their 
linearization in the small phase fluctuations yields 



c u (ip. 



s+l 



+ V's-l) 



(25) 



where c = 2 cos{2t:W / N) . In pseudomomentum repre- 
sentation — Qffc exp(z/cs), where ak satisfy the Airy 
equations 



dv? 



2u c {\ — cos k) ak 



(26) 



Since for large u the envelopes of Airy functions decay 
like ak u^^/^, then the average square of a phase step 
between nearest neighbor sites is 



^s+l 
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(27) 



Here we used Eq. and assumed zero correlation be- 
tween W and small phase fluctuations i^s+i — ips- 

The unknown constant A ~ cr^ul^'^ in Eq. (P7)) that 
measures the overall magnitude of phase fluctuations ips 
can be fixed by the initial condition a^{uc) — cr^ when 
the winding number is freezing. In this way we obtain 
scaling 



(u < N'^Uc) ^ u 



^1/2 



(28) 



soon after the freeze-out of the winding number. This 
scaling is confirmed in Fig. |4j\. On the other hand, for 
late enough times we find a saturation at 



(u > N^Uc) 



Ik 

N 



Ar2 



(29) 



which is confirmed in Fig. |6l This saturation means 
that asymptotically the phase step Og+i —Os is dominated 
by the smooth winding number 2TrW/N with negligible 
phase fiuctuations ■(/'s- 

Finally, in agreement with Ref. [l^l, the persistent 
flow (|24p in the Josephson regime is not stable for < 4 
when c = 2 cos(27rM^/A) < in Eq. (gS]). This explains 
why W'^ = for A = 4 and it is below the large- A^ linear 
fit for A^ = 8, see Fig. [3^. 




le+08 



FIG. 6. and two rescaled cr^ on = 8 sites at J = 0.1. 

For small tq the winding number W originates from a ran- 
dom walk of phase around the lattice and, consequently, 
»i N / . In contrast, for large tq there is an en- 
semble of smooth fields (f>s = exp {2niW s / N) with a ran- 
dom W and, consequently, ^ {2TTW/Ny equivalent to 



W2 f» iW7(27r)^ 



IX. 



CONCLUSION 



A ring of isolated Bose-Einstein condensates becomes 
increasingly correlated as the tunneling rate between the 



6 



condensates increases. Consequently, the initial random 
walk of phase around the ring smoothes and the vari- 
ance of its winding number decreases. However, when 
the phase becomes smooth enough, so that its small scale 
variations are no longer sufficient to let it occasionally 
hop across barriers separating different integer winding 
numbers, the winding number freezes out. As it is the 
critical smoothness that determines the trapped winding 



number, its variance does not depend on the quench time. 
Both this result and the underlying process are valid in a 
regime that allows for the "memory loss" and thus differs 
from the paradigmatic Kibble-Zurek mechanism. 
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Appendix: TWA versus exact simulation in a 
small system 

We use TWA (i.e., in effect take the limit of infinite 
average density n — > oo) before we investigate large tq. 
"Mathematically speaking", with this order of limits, 



TWA is accurate and the system gets excited from its 
instantaneous ground state for any tq. However, this 
formal statement is not satisfactory for a physicist. Af- 
ter all, we aim to investigate large but finite systems over 
large but finite timescales. So, what is of interest, is to 
determine ~ given a density (or a total number of atoms 
per site) n - what is a tq (n) when TWA breaks down. To 
this end, we have carried out a supplementary numerical 
study that, we believe, settles this issue. 



2 sites 




FIG. 7. Relative number variance at 2 sites for n — 
8,..., 1024. 

Pairs of Figures dZEl), (lOTTO)) . and (I11I12P present re- 
sults of exact simulations with n up to 1024 on 2, 3, and 
4 sites respectively. They show relative variance of an 
occupation number fig = alas, 

(30) 

as a function of tq. The same figures show results from 
TWA (black solid line) averaged over 10'* realisations. 
The exact results follow the TWA up to tq^ti) that in- 
creases with n, see Figure [T5] The five right-most data 
points on 2 sites can be fitted with 

TQ(n) = 49.5 (31) 

with error bars on the last digits. A similar fit to all five 
data points on 3 sites yields TQ{n) = 52.0 n^ '^'^. Within 



7 




the error bars these are hnear fits. Generahy, both in- 
creasing n and increasing number of sites (total number 
of particles) extend the range of validity of TWA. 

The linear fit (PT|) can be explained as follows. In the 
system of units used throughout our paper the interaction 
strength in the Hamiltonian ([ij is 

U= - . (32) 
n 

The case of 2 sites was discussed and tested by numer- 
ical simulations in Ref. Q. TWA was found exact for 
simulations times limited hy t J/U, where J was a 
constant tunneling rate. Extrapolating this result to our 
U = 1/n and time-dependent J — I/tq this inequality 
becomes 

tq < n 



and TWA should break down above Tqin) ^ n, in agree- 
ment with the linear fit ((3T|) . 
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FIG. 10. The same data as in Figure (9] Relative number 
variance at 3 sites for n — 16, 32, 64, 128, see the legend in 
Fig. El 
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FIG. 11. Relative number variance at 4 sites. 




FIG. 12. The same data as in Figure 1111 Relative number 
variance at 4 sites for n = 8, 32, see the legend in Fig. 1111 



Range of validity of TWA 
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FIG. 13. The largest range of accuracy of TWA TQ{n) as a 
function of density n for a system size of 2, 3, 4 sites. Here we 
define tq (n) as the tq when the exact variance deviates from 
the TWA by more than 20% for the first time. 



